#include "main.h"
#include "gldisplay.h" 
#include "glaxes.h"
#include "glstring.h"
#include <stdio.h>
#include <math.h>

#ifdef __APPLE__
GLint winPos_x  = 10;
GLint winPos_y  = 10;
GLint winWidth  = 920;
GLint winHeight = 600;
GLint subWidth  = 460;
GLint subHeight = 290;
GLint titHeight = 54;
#else
GLint winPos_x  = 0;
GLint winPos_y  = 0;
GLint winWidth  = 850;
GLint winHeight = 540;
GLint subWidth  = 420;
GLint subHeight = 250;
GLint titHeight = 60; 
#endif

#define NSTEP        1
#define DENSITY      0
#define VELOCITY_X   1
#define VELOCITY_Y   2
#define VELOCITY_ZX  3
#define VELOCITY_ZZ  4
#define TEMPERATURE  5
#define CLUSTER      6
#define NEW_WINDOW   7 

//coordinate system
static GLfloat xRot = 0;		                 // X-axis rotation
static GLfloat yRot = 0;  	                     // Y-axis rootation 
static GLfloat zRot = 0;                         // Z-axis rotation
static recCamera camera;
static recVec gOrigin = {0.0, 0.0, 0.0};         // coordinate origin
static GLfloat viewRange     = 26.0f;            // view range size
static GLfloat deltaAperture = 5.0f;
static GLfloat cx,cy,scale   = 1.0f;
static GLfloat shapeSize     = 50;

// various labels
GLint window_id[9];
GLint sub_display_mode[9];
int mouse_old_x, mouse_old_y, mouse_moving;
int view_model  = 1;
int point_size  = 2;
int fluid_exist = 1;
int part_exist  = 1;
int wind_number = 1;
int ntime       = 1;

GLuint ball_list;
GLuint axes_list;
/*
// lists for cluster calculation
int id[N_PT][N_PT]; // = {0};
int nf[N_PT]; 
int nv[N_PT]; 
int nc[N_PT];
*/

// light parameteres
GLfloat mat_specular[4]    = { 1.0, 1.0, 1.0, 1.0 };
GLfloat mat_diffuse[4]     = { 0.8, 0.1, 0.0, 1.0 };
GLfloat mat_shininess[]    = { 50.0 };
GLfloat light0_position[]  = { 50.0,  50.0, 50.0, 0.0 };
GLfloat light1_position[]  = { 50.0,  50.0,  50.0, 0.0 };
GLfloat light1_diffuse[]   = { 0.8, 0.1, 0.0, 1.0 };
GLfloat light1_specular[]  = { 1.0, 1.0, 1.0, 1.0 };
GLfloat light1_direction[] = { 0.0, 0.0, 1.0 };

/*
void processMenuEvents(int option) {

    int winID = glutGetWindow();

	switch (option) {
		case NEW_WINDOW : 
            //create_subWindow();
            break;
		case VELOCITY_X : 
            sub_display_mode[winID] = VELOCITY_X;
            break;
		case VELOCITY_Y : 
            sub_display_mode[winID] = VELOCITY_Y;
            break;
		case VELOCITY_ZX : 
            sub_display_mode[winID] = VELOCITY_ZX;
            break;
		case VELOCITY_ZZ : 
            sub_display_mode[winID] = VELOCITY_ZZ;
            break;
		case TEMPERATURE : 
            sub_display_mode[winID] = TEMPERATURE;
            break;
		case CLUSTER : 
            sub_display_mode[winID] = CLUSTER;
            break;
		case DENSITY : 
            sub_display_mode[winID] = DENSITY;
            break;
	}

    glutPostRedisplay();
}


void createGLUTMenus() {

	int menu;

	// create the menu and
	// tell glut that "processMenuEvents" will 
	// handle the events
	menu = glutCreateMenu(processMenuEvents);
	
	//add entries to our menu
	glutAddMenuEntry("New Window", NEW_WINDOW);
	glutAddMenuEntry("-----------------------",0);
	glutAddMenuEntry("Density", DENSITY);
	glutAddMenuEntry("Cluster ", CLUSTER);
	glutAddMenuEntry("Temperature ", TEMPERATURE);
	glutAddMenuEntry("Velocity in X on front", VELOCITY_X);
	glutAddMenuEntry("Velocity in Y on front", VELOCITY_Y);
    glutAddMenuEntry("Velocity in X on top", VELOCITY_ZX);
	glutAddMenuEntry("Velocity in Z on top", VELOCITY_ZZ);
	
	// attach the menu to the right button
	glutAttachMenu(GLUT_RIGHT_BUTTON);
	//glutAttachMenu(GLUT_LEFT_BUTTON);
}

*/
/* ------------------------------------------------------------------------*/
// functions about display
int opengl_display(int argc, char ** argv)
{
	viewRange = XSIZE / 2 + 4.0;
    shapeSize =  XSIZE;
    printf("\n %f \n", viewRange);

    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_DOUBLE|GLUT_DEPTH|GLUT_RGB);

    int ox = winPos_x;
    int oy = winPos_y;
    int ww = winWidth;
    int wh = winHeight;
    //int sw = subWidth;
    //int sh = subHeight;
    //int th = titHeight;
    create_mainWindow(ox, oy, ww, wh);
    /*

    create_subWindow(ox+ww, oy, sw, sh, DENSITY);
    create_subWindow(ox+ww, oy+sh+th, sw, sh, TEMPERATURE);
    create_subWindow(ox, oy+wh+th, sw, sh, VELOCITY_X);
    create_subWindow(ox+sw, oy+wh+th, sw, sh, VELOCITY_Y);
    create_subWindow(ox+ww, oy+wh+th, sw, sh, CLUSTER);
    */
    //PG_OneStep();
    //glutPostRedisplay();
    //glutIdleFunc(GL_Idle);
    glutMainLoop();

    return 0;
}

int create_mainWindow(int pos_x, int pos_y, int w, int h)
{

    int winID;
    char * nameOfProgram = "     MD-CPP";
    //printf("This program is for %s \n", nameOfProgram);

    glutInitWindowSize(w, h);                          // initialize main window
    glutInitWindowPosition(pos_x, pos_y);                      // specify its position
    winID = glutCreateWindow(nameOfProgram);             // create the main window  
    window_id[1] = winID;
    glutDisplayFunc(GL_Main_Display);                    // register the callbacks
    glutReshapeFunc(GL_Main_Resize);
    glutMotionFunc(GL_Motion);
    glutMouseFunc(GL_Mouse);
    glutKeyboardFunc(GL_Keyboard);
    GL_Generate_CallLists();

    return 0;
	
}

/*
int create_subWindow(int pos_x, int pos_y, int w, int h, int TYPE)
{
    int winID;

    glutInitWindowSize(w, h);                          // initialize 1st sub-window
    glutInitWindowPosition(pos_x, pos_y);                   // specify its position
    winID = glutCreateWindow("Distribution");            // create the 1st sub-window
    window_id[++wind_number] = winID; 
    sub_display_mode[winID] = TYPE;                   // specify the type of window     
    glutDisplayFunc(GL_Sub_Display);                     // register the callbacks
    glutReshapeFunc(GL_Sub_Resize);
    glutKeyboardFunc(GL_Sub_Keyboard);
    createGLUTMenus();	

    return 0;
}

*/
void GL_DrawBall()
{
    glMaterialfv(GL_FRONT, GL_DIFFUSE, mat_diffuse);
    glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular);
    glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess);
    glLightfv(GL_LIGHT0, GL_POSITION, light0_position);
    //glLightfv(GL_LIGHT1, GL_POSITION, light1_position);
    //glLightfv(GL_LIGHT1, GL_DIFFUSE, light1_diffuse);
    //glLightfv(GL_LIGHT1, GL_SPECULAR, light1_specular);
    //glLightfv(GL_LIGHT1, GL_SPOT_DIRECTION, light1_direction);
    glEnable(GL_LIGHTING);
    glEnable(GL_LIGHT0);
    //glColor4f(0.5, 0.5, 0.0, 0.8);
    glutSolidSphere(0.2, 40, 40);
	
    glDisable(GL_LIGHTING);
}

void GL_DrawBox()
{
    double x, y, z;
    double boxRatio = 1.01;	
    
	static float gray[4] = {0.6, 0.6, 0.6, 0.8};
	glColor4fv(gray);

    x = boxRatio * XSIZE;
    y = boxRatio * YSIZE *0.9;
    z = boxRatio * ZSIZE;

    glPushMatrix();
    glTranslatef((XSIZE-x)*0.5, 0.0, (ZSIZE-z)*0.5);

	
    glBegin(GL_LINE_LOOP);          //Left Surface
	glVertex3d(x, y, z);
	glVertex3d(x, 0, z);
	glVertex3d(0, 0, z);
	glVertex3d(0, y, z);
    glEnd();
	
    glBegin(GL_LINE_LOOP);          //Right Surface
	glVertex3d(x, y, 0);
	glVertex3d(x, 0, 0);
	glVertex3d(0, 0, 0);
	glVertex3d(0, y, 0);
    glEnd();
	
    glBegin(GL_LINES);
	glVertex3d(x, y, z);
	glVertex3d(x, y, 0);
	glVertex3d(x, 0, z);
	glVertex3d(x, 0, 0);
	glVertex3d(0, 0, z);
	glVertex3d(0, 0, 0);
	glVertex3d(0, y, z);
	glVertex3d(0, y, 0);
    glEnd();
    glPopMatrix();
	
    /*
	 glBegin(GL_LINES);
	 glVertex3fv(light0_position);
	 glVertex3f(0, 0, 0);
	 glEnd();
	 
	 glBegin(GL_LINES);
	 glVertex3fv(light1_position);
	 glVertex3f(0, 0, 0);
	 glEnd();
	 */
}

void GL_Init(void)
{
    //  select clearing (background) color 
    glClearColor(0.0, 0.0, 0.0, 0.0);
	
    //  set projection
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
    glOrtho(-viewRange, viewRange, -viewRange, viewRange, -viewRange, viewRange);
    glMatrixMode(GL_MODELVIEW);
	
    glEnable(GL_POINT_SMOOTH);
    glEnable(GL_BLEND);
    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
    glEnable(GL_COLOR_MATERIAL);
    glEnable(GL_DEPTH_TEST);
	
}

void GL_Main_Display(void)
{
    int i, j;
    float c;
	char str[50];
    float white[4] = {1.0, 1.0, 1.0, 1.0};
    sprintf(str, " kb = %d", 1);

	
    ///  clear all pixels 
    glDepthMask(GL_TRUE);
    glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);
    glLoadIdentity();

	// printout the  parameters
    glutPrint(0.6*viewRange, 0.8*viewRange, 0.0, str, white, GLUT_BITMAP_TIMES_ROMAN_24);
    
	// draw the axis 
	glPushMatrix();
	glTranslatef(-0.92*viewRange, -0.88*viewRange, 0);
	glRotatef(yRot, 1, 0, 0);
	glRotatef(xRot, 0, 1, 0);
	glRotatef(zRot, 0, 0, 1);
	glScalef(scale,scale*winWidth/winHeight,scale);
	glCallList(axes_list);
    glutPrint(3.5f, -0.2f, 0.0f, "X", white, GLUT_BITMAP_TIMES_ROMAN_10);	
    glutPrint(-0.15f, 3.4f, 0.0f, "Y", white, GLUT_BITMAP_TIMES_ROMAN_10);	
    glutPrint(0.6f, -0.6f, 3.2f, "Z", white, GLUT_BITMAP_TIMES_ROMAN_10);	
    glPopMatrix();
	
	// draw the main screnio
    glRotatef(yRot, 1, 0, 0);        
    glRotatef(xRot, 0, 1, 0);
    glRotatef(zRot, 0, 0, 1);
	
	//draw the box
    glPushMatrix();
	glTranslatef(-0.5*shapeSize*scale+cx, -0.5*shapeSize*scale+cy+4, -0.5*shapeSize*scale);
	glScalef(scale,scale*winWidth/winHeight,scale);
	GL_DrawBox();

	//draw the particles
	glPointSize(point_size);


	if(fluid_exist)
		j = 0;
	else
		j = 1;
	for (; j < NS; j++)  {
		if(part_exist && j == 2) continue;
		c = ( 0.9 - (j+1) * 0.4 ) * 255;
		glColor3ub((GLubyte) c, (GLubyte) (128+c*0.5), 255);
		if(j == 2) glColor3ub(255, 25, 24);
		glBegin(GL_POINTS);
		for(i = LL[j]-1; i < LU[j]; i++) {
		    //if( i%10==9) glColor3ub(255, 25, 24);
		    //if( i==23||i==24) glColor3ub(255, 25, 24);
            //else glColor3ub((GLubyte) c, (GLubyte) (128+c*0.5), 255);
			glVertex3f((*R0[i]).x,(*R0[i]).y, (*R0[i]).z);      
		}
		glEnd();
	}
	
    glPopMatrix();

    /*
	for (; j < NS-1; j++)  {
		for (i = LL[j]-1; i < lu[j]; i++) {
			c = ( 0.9 - myParticles[i].id * 0.4 ) * 255;
			glColor3ub((GLubyte) c, (GLubyte) (128+c*0.5), 255);        
			//if(myParticles[i].ic == 4) glColor3ub(255, 150, 0);
			//if(myParticles[i].ic == 5) glColor3ub(150, 25, 240);
			//if(myParticles[i].ic == 5) glColor3ub(250, 25, 24);
			if(myParticles[i].id == 4) glColor3ub(25, 25, 245);
			glPushMatrix();
			glTranslatef(myParticles[i].position[0], myParticles[i].position[1], myParticles[i].position[2]);
			glCallList(ball_list);
			glPopMatrix();
		}
	}
	
	 glPushMatrix();
	 glPointSize(point_size);
	 glTranslatef(-0.5*shapeSize*scale+cx, -0.5*shapeSize*scale+cy+3, -0.5*shapeSize*scale);
	 glScalef(scale,scale*winWidth/winHeight,scale);
	 if(fluid_exist == 0)
	 j = 0;
	 else
	 j = 1;
	 for (; j < NS-3; j++)  {
	 if(part_exist == 0 && j == 2) continue;
	 c = ( 0.9 - (j+1) * 0.4 ) * 255;
	 glColor3ub((GLubyte) c, (GLubyte) (128+c*0.5), 255);
	 if(j == 2) glColor3ub(25, 250, 240);
	 glBegin(GL_POINTS);
	 for (i = LL[j]-1; i < lu[j]; i++) {
	 //c = ( 0.9 - myParticles[i].ic * 0.4 ) * 255;
	 //glColor3ub((GLubyte) c, (GLubyte) (128+c*0.5), 255);
	 //if(myParticles[i].ic == 3) glColor3ub(25, 250, 240);
	 glVertex3fv(myParticles[i].position);      
	 }
	 glEnd();
	 }
	 
	 for (; j < NS-1; j++)  {
	 for (i = LL[j]-1; i < lu[j]; i++) {
	 c = ( 0.9 - myParticles[i].id * 0.4 ) * 255;
	 glColor3ub((GLubyte) c, (GLubyte) (128+c*0.5), 255);        
	 //if(myParticles[i].ic == 4) glColor3ub(255, 150, 0);
	 //if(myParticles[i].ic == 5) glColor3ub(150, 25, 240);
	 //if(myParticles[i].ic == 5) glColor3ub(250, 25, 24);
	 if(myParticles[i].id == 4) glColor3ub(25, 25, 245);
	 glPushMatrix();
	 glTranslatef(myParticles[i].position[0], myParticles[i].position[1], myParticles[i].position[2]);
	 glCallList(ball_list);
	 glPopMatrix();
	 }
	 }
	 glPopMatrix();
	*/

    glutSwapBuffers();
}

/*
void GL_Sub_Display(void)
{
	int i;
    int winID = glutGetWindow();
    int winDM = sub_display_mode[winID];
    //printf("current window is %d\n", winID);
	static int imax = 40;
    static float gray[4] = {0.6, 0.6, 0.6, 0.8};
    static float rose[4] = {0.9, 0.2, 0.5, 0.8};
    float  distY[NZ];
	GLfloat bar;
    GLfloat axL = 0.92;
    GLfloat trL = 0.03;
    char str[8];


	glClear(GL_COLOR_BUFFER_BIT);               // clear all pixels 
    
	glScalef(1.0,subWidth/subHeight,1.0);
    glLineWidth(1.0);


    if(winDM == CLUSTER)                        // graph for cluster
    {
    glutSetWindowTitle("Cluster Distribution");
    bar =  axL / imax * 0.95;

    glPushMatrix();
    glTranslatef(0.02, 0.03, 0);                // origint setup
    glBegin(GL_LINES);
		//glColor3f(0.2, 0.2, 0.4);		
        glColor4fv(gray);
        glVertex3f(0.0, 0.0, 0.0);
        glVertex3f(0.0, axL, 0.0);
        glVertex3f(axL, 0.0, 0.0);
        glVertex3f(0.0, 0.0, 0.0);
    glEnd();
    
    glBegin(GL_POLYGON);
		glVertex3f(trL/3, axL-trL, 0.0);
		glVertex3f(-trL/3, axL-trL, 0.0);
		glVertex3f(0.0, axL, 0.0);
	glEnd();

    glBegin(GL_POLYGON);
		glVertex3f(axL-trL, trL/2, 0.0);
		glVertex3f(axL-trL, -trL/2, 0.0);
		glVertex3f(axL, 0.0, 0.0);
	glEnd();

    for(i = 1; i <= imax; i++)
    {   
        glBegin(GL_LINES);
        glVertex3f((i-0.5)*bar, 0.0, 0.0);
        glVertex3f((i-0.5)*bar, 0.02*axL, 0.0);
        glEnd();
        if(i%5==0)
        {
            sprintf(str, "%d", i);
            glutPrint((i-0.8)*bar, -0.04, 0.0, str, gray, GLUT_BITMAP_TIMES_ROMAN_10);
        }
    }

    glColor4fv(rose);
    for(i = 1; i <= imax; i++)
    {    
        glBegin(GL_POLYGON);
        glVertex3f((i-1)*bar, 0.0, 0.0);
        glVertex3f((i-1)*bar, nc[i]*0.1, 0.0);
        glVertex3f(i*bar, nc[i]*0.1, 0.0);
        glVertex3f(i*bar, 0.0, 0.0);
        glEnd();
    }

    glPopMatrix();
    }
    else
    {
           
    int i, j; 
    int nx = NX;
    int ny = NY;
    int nz = NZ;
    float ratio;
    GLfloat inx = 1.0 / nx;
    GLfloat iny = 1.0 / ny;

    if(winDM==DENSITY)
    {
        load_data("./Density.dis", nx, ny);
        glutSetWindowTitle("Density Distribution");
        if(dataMax > 0.0065) dataMax = 0.0065;
    }
    if(winDM==TEMPERATURE)
    {
        load_data("./Temp.dis", nx, ny);
        glutSetWindowTitle("Temperature Distribution");
        if(dataMax > 2.0) dataMax = 2.0;
    }
    if(winDM==VELOCITY_X) 
    {
        load_data("./VeloX.dis", nx, ny);
        glutSetWindowTitle("Velocity in X on front");
        if( dataMax > 2.0) dataMax = 2.0;
        if( dataMax > -dataMin ) dataMin = -dataMax;
        dataMax = -dataMin;
    }
    if(winDM==VELOCITY_Y)
    {
        load_data("./VeloY.dis", nx, ny);
        glutSetWindowTitle("Velocity in Y on front");
        if( dataMax > 2.0) dataMax = 2.0;
        if( dataMax > -dataMin ) dataMin = -dataMax;
        dataMax = -dataMin;
    }
    if(winDM==VELOCITY_ZX)
    {
        ny = NZ;
        iny = 1.0 / ny;
        load_data("./VelZX.dis", nx, ny);
        glutSetWindowTitle("Velocity in X on top");
        if( dataMax > 2.0) dataMax = 2.0;
        if( dataMax > -dataMin ) dataMin = -dataMax;
        dataMax = -dataMin;
    }
    if(winDM==VELOCITY_ZZ)
    {
        ny = NZ;
        iny = 1.0 / ny;
        load_data("./VelZZ.dis", nx, nz);
        glutSetWindowTitle("Velocity in Z on top");
        if( dataMax > 2.0) dataMax = 2.0;
        if( dataMax > -dataMin ) dataMin = -dataMax;
        dataMax = -dataMin;
    }



    glColor4f(0.5, 0.5, 0.5, 0.8);             //draw grids lines
    glBegin(GL_LINES);
        for(i = 0; i <= nx; i++ )
        {
            glVertex3f(inx * i, 0.0, 0.0);
            glVertex3f(inx * i, 1.0, 0.0);
        }
        for(j = 0; j <= ny; j++ )
        {
            glVertex3f(0.0, iny * j, 0.0);
            glVertex3f(1.0, iny * j, 0.0);
        }
    glEnd();



    glPushMatrix();
    glTranslatef(-inx, 1.0, 0.0);
    for(i =0; i < nx; i++)
    {
        glTranslatef(inx, -1.0, 0.0);
        for(j = 0; j < ny; j++)
        {
            ratio = (dist[j][i]) / (dataMax);
            glColor4f(0.4*ratio, 0.4*ratio, 1.0*ratio, 0.9);
            if(dataMin < 0.0)
            {
                glColor4f(0.2*ratio, 1.0*ratio, 0.4*ratio, 0.9);
            }
            if(ratio < 0.0)
            {
                ratio = (dist[j][i] ) / ( dataMin );
                glColor4f(1.0*ratio, 0.2*ratio, 0.0*ratio, 0.9);
            }
            glBegin(GL_POLYGON);
                glVertex3f(0.0, 0.0, 0.0);
                glVertex3f(inx, 0.0, 0.0);
                glVertex3f(inx, iny, 0.0);
                glVertex3f(0.0, iny, 0.0);
            glEnd();
            glTranslatef(0.0, iny, 0.0);
        }
    }
    glPopMatrix();


        dataMax = 0.0;
        dataMin = 0.0;
        for(j = 0; j < ny; j++)
        {
           distY[j] = 0;
           for(i =0; i < nx; i++)
           { 
                distY[j]+=dist[j][i];
           } 
           distY[j]/=nx;
           if(distY[j]>dataMax)dataMax = distY[j];
           if(distY[j]<dataMin)dataMin = distY[j];

        }

        glPushMatrix();                               
        glTranslatef(1.0+inx, -iny, 0.0);
        for(j = 0; j < ny; j++)
        {
           glTranslatef(0.0, iny, 0.0);
           ratio = (distY[j]) / (dataMax);
           glColor4f(0.4*ratio, 0.4*ratio, 1.0*ratio, 0.9);
               if(dataMin < 0.0)
               { 

               if(distY[j] > 0.0)
               { 
                    ratio = (distY[j]) / (dataMax);
                    glColor4f(0.2*ratio, 1.0*ratio, 0.4*ratio, 0.9);
               } 
               if(distY[j] < 0.0)
               { 
                    ratio = (distY[j] ) / ( dataMin );
                    glColor4f(1.0*ratio, 0.2*ratio, 0.0*ratio, 0.9);
               } 

               }
               glBegin(GL_POLYGON);
                    glVertex3f(0.0, 0.0, 0.0);
                    glVertex3f(inx, 0.0, 0.0);
                    glVertex3f(inx, iny, 0.0);
                    glVertex3f(0.0, iny, 0.0);
               glEnd();
        }
        glPopMatrix();

        dataMax = 0.0;
        dataMin = 0.0;
        for(i = 0; i < nx; i++)
        {
           distY[i] = 0;
           for(j = 0; j < ny; j++)
           { 
                distY[i]+=dist[j][i];
           } 
           distY[i]/=ny;
           if(distY[i]>dataMax)dataMax = distY[i];
           if(distY[i]<dataMin)dataMin = distY[i];

        }

        glPushMatrix();                               
        glTranslatef( 0.0, -iny*1.2, 0.0);
        for(i = 0; i < nx; i++)
        {
           glTranslatef(inx, 0.0, 0.0);
           ratio = (distY[i]) / (dataMax);
           glColor4f(0.4*ratio, 0.4*ratio, 1.0*ratio, 0.9);
               if(dataMin < 0.0)
               { 

               if(distY[i] > 0.0)
               { 
                    ratio = (distY[i]) / (dataMax*1.4);
                    glColor4f(0.2*ratio, 1.0*ratio, 0.4*ratio, 0.9);
               } 
               if(distY[i] < 0.0)
               { 
                    ratio = (distY[i] ) / ( dataMin*1.4 );
                    glColor4f(1.0*ratio, 0.2*ratio, 0.0*ratio, 0.9);
               } 

               }
               glBegin(GL_POLYGON);
                    glVertex3f(0.0, 0.0, 0.0);
                    glVertex3f(inx, 0.0, 0.0);
                    glVertex3f(inx, iny, 0.0);
                    glVertex3f(0.0, iny, 0.0);
               glEnd();
        }
        glPopMatrix();

    glLineWidth(3.0);
    if(winDM==VELOCITY_X)                        // draw for legend
    {
        axL = 0.06;
        glPushMatrix();
            glTranslatef(0.12, 0.9, 0);
            glBegin(GL_LINES);
                glColor4f(0.0, 0.0, 0.0, 0.9);
                glVertex3f(0.0, 0.0, 0.0);
                glColor4f(0.2, 1.0, 0.4, 0.9);
                glVertex3f(axL, 0.0, 0.0);
            glEnd();
        
            glBegin(GL_POLYGON);
	        	glVertex3f(axL, +trL/2, 0.0);
	        	glVertex3f(axL, -trL/2, 0.0);
	        	glVertex3f(axL+trL, 0.0, 0.0);
       	    glEnd();

            glBegin(GL_LINES);
                glColor4f(0.0, 0.0, 0.0, 0.9);
                glVertex3f(0.0, 0.0, 0.0);
                glColor4f(1.0, 0.2, 0.0, 0.9);
                glVertex3f(-axL, 0.0, 0.0);
            glEnd();


            glBegin(GL_POLYGON);
	        	glVertex3f(-axL, +trL/2, 0.0);
	        	glVertex3f(-axL, -trL/2, 0.0);
	        	glVertex3f(-axL-trL, 0.0, 0.0);
       	    glEnd();
        glPopMatrix();
    }

    if(winDM==VELOCITY_ZX)                     // draw for legend
    {
        axL = 0.06;
        glPushMatrix();
            glTranslatef(0.12, 0.9, 0);
            glBegin(GL_LINES);
                glColor4f(0.0, 0.0, 0.0, 0.9);
                glVertex3f(0.0, 0.0, 0.0);
                glColor4f(0.2, 1.0, 0.4, 0.9);
                glVertex3f(axL, 0.0, 0.0);
            glEnd();
        
            glBegin(GL_POLYGON);
	        	glVertex3f(axL, +trL/2, 0.0);
	        	glVertex3f(axL, -trL/2, 0.0);
	        	glVertex3f(axL+trL, 0.0, 0.0);
       	    glEnd();

            glBegin(GL_LINES);
                glColor4f(0.0, 0.0, 0.0, 0.9);
                glVertex3f(0.0, 0.0, 0.0);
                glColor4f(1.0, 0.2, 0.0, 0.9);
                glVertex3f(-axL, 0.0, 0.0);
            glEnd();


            glBegin(GL_POLYGON);
	        	glVertex3f(-axL, +trL/2, 0.0);
	        	glVertex3f(-axL, -trL/2, 0.0);
	        	glVertex3f(-axL-trL, 0.0, 0.0);
       	    glEnd();
        glPopMatrix();
    }


    if(winDM==VELOCITY_Y) 
    {
        axL = 0.08;
        glPushMatrix();
            glTranslatef(0.05, 0.84, 0);
            glBegin(GL_LINES);
                glColor4f(0.0, 0.0, 0.0, 0.9);
                glVertex3f(0.0, 0.0, 0.0);
                glColor4f(0.2, 1.0, 0.4, 0.9);
                glVertex3f(0.0, axL, 0.0);
            glEnd();
        
            glBegin(GL_POLYGON);
	        	glVertex3f( trL/3, axL, 0.0);
	        	glVertex3f(-trL/3, axL, 0.0);
	        	glVertex3f(0.0, axL+trL, 0.0);
       	    glEnd();

            glBegin(GL_LINES);
                glColor4f(0.0, 0.0, 0.0, 0.9);
                glVertex3f(0.0, 0.0, 0.0);
                glColor4f(1.0, 0.2, 0.0, 0.9);
                glVertex3f(0.0, -axL, 0.0);
            glEnd();


            glBegin(GL_POLYGON);
	        	glVertex3f( trL/3, -axL, 0.0);
	        	glVertex3f(-trL/3, -axL, 0.0);
	        	glVertex3f(0.0, -axL-trL, 0.0);
       	    glEnd();
        glPopMatrix();
    }

    if(winDM==VELOCITY_ZZ) 
    {
        axL = 0.08;
        glPushMatrix();
            glTranslatef(0.05, 0.84, 0);
            glBegin(GL_LINES);
                glColor4f(0.0, 0.0, 0.0, 0.9);
                glVertex3f(0.0, 0.0, 0.0);
                glColor4f(0.2, 1.0, 0.4, 0.9);
                glVertex3f(0.0, axL, 0.0);
            glEnd();
        
            glBegin(GL_POLYGON);
	        	glVertex3f( trL/3, axL, 0.0);
	        	glVertex3f(-trL/3, axL, 0.0);
	        	glVertex3f(0.0, axL+trL, 0.0);
       	    glEnd();

            glBegin(GL_LINES);
                glColor4f(0.0, 0.0, 0.0, 0.9);
                glVertex3f(0.0, 0.0, 0.0);
                glColor4f(1.0, 0.2, 0.0, 0.9);
                glVertex3f(0.0, -axL, 0.0);
            glEnd();


            glBegin(GL_POLYGON);
	        	glVertex3f( trL/3, -axL, 0.0);
	        	glVertex3f(-trL/3, -axL, 0.0);
	        	glVertex3f(0.0, -axL-trL, 0.0);
       	    glEnd();
        glPopMatrix();
    }

    }

    glutSwapBuffers();

}
*/
void GL_Sub_Resize(int width, int height)
{
    glViewport(0, 0, width, height);
/*	select clearing (background) color */
	glClearColor(0.0, 0.0, 0.0, 0.0);

/*	initialize viewing values */
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
//	glOrtho(0.0, 1.0, 0.0, 1.0, -1.0, 1.0);
	gluOrtho2D(-0.1, 1.1, -0.1, 1.1);

    glEnable(GL_LINE_SMOOTH);
    glEnable(GL_POLYGON_SMOOTH);
    glEnable(GL_BLEND);
    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

	//glShadeModel(GL_FLAT);
}

void GL_Main_Resize(int width, int height)
{
	
	GLdouble ratio, radians, wd2;
	GLdouble left, right, top, bottom, zNear, zFar;
	
    glViewport(0, 0, width, height);
	
    //  set camera 
    GL_SetCam(&camera);
    camera.viewHeight = height;
    camera.viewWidth = width;
	
    // set projection
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	
	if(view_model == 1) {
	    glOrtho(-viewRange, viewRange, -viewRange, viewRange, -1.8*viewRange, 1.8*viewRange);
	}
	else {
        glOrtho(-viewRange, viewRange, -0.9*viewRange, 0.8*viewRange, -viewRange, viewRange);
	    zNear = - camera.viewPos.z - shapeSize * 0.5;
	    if (zNear < 0.00001) zNear = 0.00001;
	    zFar = - camera.viewPos.z + shapeSize * 0.5;
        if (zFar < 1.0) zFar = 1.0;
		
        radians = 0.0174532925 * camera.aperture / 2; // half aperture degrees to radians 
        wd2 = zNear * tan(radians);
        ratio = camera.viewWidth / (float) camera.viewHeight;
	    if (ratio >= 1.0) {
		    left  = -ratio * wd2;
		    right = ratio * wd2;
		    top = wd2;
		    bottom = -wd2;	
		} else {
		    left  = -wd2;
		    right = wd2;
		    top = wd2 / ratio;
		    bottom = -wd2 / ratio;	
		}
		
		glFrustum(left, right, bottom, top, zNear, zFar);
		
    }
    glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	
	
    glShadeModel(GL_SMOOTH);
    glEnable(GL_POINT_SMOOTH);
    glEnable(GL_LINE_SMOOTH);
    glEnable(GL_BLEND);
    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
    glEnable(GL_COLOR_MATERIAL);
    glEnable(GL_DEPTH_TEST);
    glEnable(GL_NORMALIZE);
	
}


/* ------------------------------------------------------------------------*/
// functions about control: Keyboard, Mouse, Camera, ...
void GL_Mouse(int button, int state, int x, int y)
{
    if( state == GLUT_DOWN) 
    {
        mouse_old_x = x;
        mouse_old_y = y;
        mouse_moving = 1;
    }
	
    if( state == GLUT_UP)
    { 
        mouse_moving = 0;
    }
}

void GL_Motion(int x, int y)
{
    if(mouse_moving)
    {
        xRot += 0.2*(x - mouse_old_x);
        yRot += 0.2*(y - mouse_old_y);
		
        mouse_old_x = x;
        mouse_old_y = y;
		
        glutPostRedisplay();
    }
}

void GL_Keyboard(unsigned char key, int x, int y) 
{
    int i;
    static int fullscreen = 0;
    static int old_x;
    static int old_y;
    static int old_width;
    static int old_height;
	
    //printf("This is Keyboard(%c)\n",key);
	
    switch (key) {
		case ESC:
			PG_Bail(0);
			break;
			
		case 'P':
			point_size++;
			glPointSize(point_size);
			break;
			
		case 'p':
			point_size--;
			if (point_size < 1)
				point_size = 1;
			glPointSize(point_size);
			break;
			
		case 'r':
			xRot = 0;
			yRot = 0;
			cx = 0;
			cy = 0;
			scale  = 1.0;
			break;
			
		case 'v':
			view_model = 1;
			old_width = glutGet(GLUT_WINDOW_WIDTH);
			old_height = glutGet(GLUT_WINDOW_HEIGHT);
			GL_Main_Resize(old_width, old_height);
			break;
			
		case 'V':
			view_model = 2;
			old_width = glutGet(GLUT_WINDOW_WIDTH);
			old_height = glutGet(GLUT_WINDOW_HEIGHT);
			GL_Main_Resize(old_width, old_height);
			break;
			
		case 'z':
			camera.aperture -= deltaAperture;
			if (camera.aperture < 0.1) // do not let aperture <= 0.1
				camera.aperture = 0.1;
			old_width = glutGet(GLUT_WINDOW_WIDTH);
			old_height = glutGet(GLUT_WINDOW_HEIGHT);
			GL_Main_Resize(old_width, old_height);
			break;
			
		case 'Z':
			camera.aperture -= deltaAperture;
			if (camera.aperture > 179.9) // do not let aperture >= 180
				camera.aperture = 179.9;
			old_width = glutGet(GLUT_WINDOW_WIDTH);
			old_height = glutGet(GLUT_WINDOW_HEIGHT);
			GL_Main_Resize(old_width, old_height);
			break;
			
		case 'q':
			scale += 0.1;
			break;
			
		case 'e':
			scale -= 0.1;
			if(scale<0.1)scale = 0.1;
			break;
			
		case 'f':
			if(fluid_exist == 0)
				fluid_exist = 1;
			else
				fluid_exist = 0;
			break;
			
		case 'F':
			if(part_exist == 0)
				part_exist = 1;
			else
				part_exist = 0;
			break;
			
		case 'b':
			ntime--;
			if(ntime < 0)ntime = 0;
            PG_OneStep();
			break;
			
		case 'n':
			ntime+=NSTEP;
            PG_OneStep();
			break;
			
		case 'w':									// Is W Key Being Pressed?
			cy+=0.6f;								// Move Object Up
			break;
			
		case 's':									// Is S Key Being Pressed?
			cy-=0.6f;								// Move Object Down
			break;
			
		case 'd':									// Is D Key Being Pressed?
			cx+=0.6f;								// Move Object Right
			break;
			
		case 'a':									// Is A Key Being Pressed?
			cx-=0.6f;								// Move Object Left
			break;
/*
        case 'c':
            create_subWindow(1500, winPos_y, 400, 360, DENSITY);
            break;
*/			
		case '~':
			fullscreen = !fullscreen;
			if (fullscreen) {
				old_x = glutGet(GLUT_WINDOW_X);
				old_y = glutGet(GLUT_WINDOW_Y);
				old_width = glutGet(GLUT_WINDOW_WIDTH);
				old_height = glutGet(GLUT_WINDOW_HEIGHT);
				glutFullScreen();
			} else {
				glutReshapeWindow(old_width, old_height);
				glutPositionWindow(old_x, old_y);
			}
			break;
    }

    for(i = 1; i <= wind_number; i++)
    {
            glutPostWindowRedisplay(window_id[i]);
    }
}

/*
void GL_Sub_Keyboard(unsigned char key, int x, int y) 
{
    int i;
    //GLint currentWindow;
	
    switch (key) {
		case ESC:
            //currentWindow = glutGetWindow();
            //glutDestroyWindow(currentWindow); 
			PG_Bail(0);
			break;
					
		case 'P':
			point_size++;
			glPointSize(point_size);
			break;
			
		case 'p':
			point_size--;
			if (point_size < 1)
				point_size = 1;
			glPointSize(point_size);
			break;
			
		case 'r':
			xRot = 0;
			yRot = 0;
			cx = 0;
			cy = 0;
			scale  = 1.0;
			break;
			
		case 'q':
			scale += 0.1;
			break;
			
		case 'e':
			scale -= 0.1;
			if(scale<0.1)scale = 0.1;
			break;
			
		case 'f':
			if(fluid_exist == 0)
				fluid_exist = 1;
			else
				fluid_exist = 0;
			break;
			
		case 'F':
			if(part_exist == 0)
				part_exist = 1;
			else
				part_exist = 0;
			break;
			
        case 'b':
			ntime--;
			if(ntime < 0)ntime = 0;
            PG_OneStep();
			break;
			
		case 'n':
			ntime+=NSTEP;
            PG_OneStep();
			break;
            			
		case 'w':									// Is W Key Being Pressed?
			cy+=0.6f;								// Move Object Up
			break;
			
		case 's':									// Is S Key Being Pressed?
			cy-=0.6f;								// Move Object Down
			break;
			
		case 'd':									// Is D Key Being Pressed?
			cx+=0.6f;								// Move Object Right
			break;
			
		case 'a':									// Is A Key Being Pressed?
			cx-=0.6f;								// Move Object Left
			break;

        case 'c':
            create_subWindow(winPos_x, winPos_y, subWidth, subHeight, DENSITY);
            break;
	
    }
 
    for(i = 1; i <= wind_number; i++)
    {
            glutPostWindowRedisplay(window_id[i]);
    }
}
*/
void GL_SetCam(recCamera *cam)
{
    cam->aperture = 40;
    cam->rotPoint = gOrigin;
	
    cam->viewPos.x = 0.0;
    cam->viewPos.y = 0.0;
    cam->viewPos.z = -40.0;
    cam->viewDir.x = -cam->viewPos.x; 
    cam->viewDir.y = -cam->viewPos.y; 
    cam->viewDir.z = -cam->viewPos.z;
	
    cam->viewUp.x = 0;  
    cam->viewUp.y = 1; 
    cam->viewUp.z = 0;
	
    cam->viewHeight = winHeight;
    cam->viewWidth = winWidth;
}


/* ------------------------------------------------------------------------*/
// program procedure functions
void PG_Bail(int code)
{
    int i;

    glDeleteLists(ball_list, 1);
    //free(myParticles);
    //printf("This is bail(%d)\n",code);
    //printf("The Program is going to end~\n");
    /* shut down our window */
    for(i = 1; i <= wind_number; i++)
    {
            glutDestroyWindow(window_id[i]);
    }
    /* exit the program...normal termination. */
    exit(code);
}
/*
void PG_Cluster(psParticle* p)
{
    int i, j,  ipar, jpar, find;
    int idx, idy, idz;
    int lpar = 32;
	int npar = (lu[4]-lu[2])/(lpar);
    real dx, dy, dz, xi2, yi2, zi2;
    real dd, d2 = 1;                   //Cutoff
	
    xi2 = 2.0 / XSIZE;
    yi2 = 2.0 / YSIZE;
    zi2 = 2.0 / ZSIZE;
    
    for(i=0;i<N_PT;i++)
    {
        nf[i] = 0;
        nv[i] = 0;
        nc[i] = 0;
        for(j=0;j<N_PT;j++)
            id[i][j] = 0;
    }
	// for(i=0;i<10;i++)
	// for(j=0;j<10;j++)
	// printf("%d", id[i][j]);
	 
	
    //printf("\tnpar = %d lpar = %d total = %d\n", npar, lpar, lu[4]-lu[2]);
    if(npar*lpar != lu[4]-lu[2])
    {
        printf("bad particle count!\n");
        exit(0);
    }
    
	
    for(ipar = 0; ipar < npar - 1; ipar++)
    {
        for(jpar = ipar+1; jpar < npar; jpar++)
        {
			
            find = 0;
            if(find == 1) break;
            for(i = lu[2] + ipar * lpar; i < lu[2] + (ipar+1) * lpar; i++)
            {
                if(find == 1) break;
                for(j = lu[2] + jpar * lpar; j < lu[2] + (jpar+1) * lpar; j++)
                {
                    dx = p[i].position[0] - p[j].position[0];
                    dy = p[i].position[1] - p[j].position[1];
                    dz = p[i].position[2] - p[j].position[2];
                    idx = dx * xi2;
                    idy = dy * yi2;
                    idz = dz * zi2;
                    dx -= idx * XSIZE;
                    dy -= idy * YSIZE;
                    dz -= idz * ZSIZE;
                    dd = dx * dx + dy * dy + dz *dz;
                    if( dd < d2)
                    {
                        id[jpar+1][ipar+1] = 1;
                        //printf("id[%d][%d] = %f\n", jpar+1, ipar+1, dd);
                        find = 1;
                        break;
                    }
                }
            }
        }
    }
    
    find = 0;
    for(i=0;i<npar;i++)
        for(j=0;j<npar;j++)
            find += id[i][j];
    //printf("%d pair of myParticles found\n", find);
	
	
    nf[1] = 1;
    for(i = 2; i <= npar; i++)
    {
        nf[i] = i;
        for(j = 1; j <= i - 1; j++)
        {
            nf[j] = nf[nf[j]];
            if(id[i][j] == 1) nf[nf[nf[j]]] = i;
        }
		
    }
	
    for(i = 1; i <= npar; i++)
    {
        nf[i] = nf[nf[i]];
    }
	
    for(i = 1; i <= npar; i++)
    {
        j = nf[i];
        nv[j]++;
    }
	
    for(i = 1; i <= npar; i++)
    {
        j = nv[i];
        if(j > 0) nc[j]++;
    }
	
    for(find = npar; npar > 0; find--)
        if(nc[find] > 0) break;
	
    printf("\t[%d]\t",ntime);
    for(i = 1; i <= find; i++)
        printf("%d  ", nc[i]);
	
    printf("\n");
}
*/
void PG_OneStep()
{
	//load_particles("./fort.22", myParticles);
	//PG_Cluster(myParticles);
	glutPostRedisplay();
	
}

void GL_Idle()
{
	if(ntime < 1000)
	{
		ntime++;
		//load_particles("./fort.22", myParticles);
		//PG_Cluster(myParticles);
		glutPostRedisplay();
	}
	else
	{
		exit(0);
	}
	
}

void GL_Generate_CallLists()
{
    axes_list = glGenLists(1);
    ball_list = glGenLists(2);

    glNewList(axes_list, GL_COMPILE);
    GL_DrawAxes(viewRange);
    glEndList();

    glNewList(ball_list, GL_COMPILE);
    GL_DrawBall();
    glEndList();
}

